Temporal changes in zooplankton indicators highlight a bottom-up process in the Bay of Marseille (NW Mediterranean Sea)

Sixteen years (2005–2020) of zooplankton monitoring in the Bay of Marseille (N-W Mediterranean Sea) are analyzed in relation to physical, meteorological, climatic and biotic data. Samples were collected every two weeks by a vertical haul (0–55 m) of a 200 μm plankton net. Different indices characterizing the mesozooplankton are compared: biomass dry weight of four size fractions between 200 and 2000 μm; abundances of the whole of the mesozooplankton and of 13 main taxonomic groups defined from plankton imagery; seasonal onset timing of each zooplankton group; and two other types of indices: the first characterized diversity based on abundance data, and the second was derived from zooplankton size spectra shape. The clearest pattern in the environmental compartment was an overall decreasing trend in nutrients, shifts in phytoplankton metrics (i.e. size structure and particulate organic matter), and changes in winter conditions (i.e. increasing temperatures, precipitation and NAO). Interannual patterns in the mesozooplankton community were: (i) a decrease of total abundance (ii) a decrease in biomass for the four size fractions, with an earlier decrease for the 1000–2000 μm size fraction (in 2008); (iii) a reduced dominance of copepods (calanoids and oithonoids) and a concomitant increase in abundance of other taxonomic groups (crustaceans, pteropods, chaetognaths, salps) which induced higher diversity; (iv) a first shift in size spectra towards smaller sizes in 2009, when the 1000–2000 μm size fraction biomass decreased, and a second shift towards larger sizes in 2013 along with increased diversity; and (iv) a later onset in the phenology for some zooplankton variables and earlier onset for salps. Concomitant changes in the phytoplankton compartment, winter environmental conditions, zooplankton community structure (in size and diversity) and zooplankton phenology marked by a shift in 2013 suggest bottom-up control of the pelagic ecosystem.


1-Introduction
Zooplankton are short-lived organisms, highly sensitive to environmental changes, which makes their long-term monitoring an excellent approach to provide responses to management issues [1,2].Previous time series analyses have shown that zooplankton shifts can be observed at various scales: at large scale through climate-related processes [3,4], at regional scale in relation with regional meteorological phenomena, and at local scale (e.g.bays) most often in response to human pressure [5,6].These multiple stressors alter the physiology, behavior and phenology of organisms, and consequently population distributions [7,8].Nonetheless the multiple scales of forcing make it difficult to disentangle the causation.
The increase in zooplankton observations has made it possible to validate several indicators such as those based on individual taxa, species assemblages, diversity indices or size spectra, which are useful for characterizing the state of pelagic ecosystems [1].The development of multi-indicator approaches has made it possible to better compare them, and thus to identify their complementary contributions.Therefore nowadays quantification of the roles of zooplankton in biogeochemical fluxes and trophic fluxes towards higher levels requires, in addition to information on biomass and size spectra, knowledge of taxonomic distributions and associated ecological traits.
In the Gulf of Lion (NW Mediterranean), earlier studies dedicated to zooplankton were focused on quantification of either biogeochemical flux or pelagic trophic flows, and were all carried out during oceanographic cruises over short periods (less than a month) [9][10][11][12][13][14][15][16][17][18][19].Recent increased interest in the role of zooplankton as trophic resources for planktivorous fish has been prompted by the reduction in size of two planktivorous fishes (European pilchard, Sardina pilchardus, and European anchovy, Engraulis encrasicolus) in the Gulf of Lion in 2008 (NW Mediterranean) [20][21][22][23].Various studies have led to rejection of the hypotheses of the impact of overfishing, increased predation, disease and a decline in recruitment [24].The observation of a drastic reduction in the size of individuals [20,24,25] has been linked to a deficit in individual growth attributed to potential changes in plankton resources, although this has not yet been proven.To shed light on this potential planktonic change, the data from oceanographic cruises over the last few decades in the Gulf of Lion are too few and limited in time, and the closest locations of zooplankton monitoring time series in the Ligurian Sea [4,26,27] or the Balearic Sea [28] are far apart.
A mesozooplankton monitoring survey that has been carried out since 2005 in the Bay of Marseille (BoM) provides a zooplankton time series in the Gulf of Lion which encompasses the critical shift period.Zooplankton sampling was added to an existing monitoring programme started in 1994 including hydrological and nutrient measurements, as well as sampling for a few biotic variables (pigment, bacteria and phytoplankton).The sampling station at the center of the BoM is subject to multiple regional environmental influences such as winter deep convection phenomena from the Provencal basin off the Gulf of Lion, north-westerly winds (known as Mistral), which induce coastal upwelling events [29], intrusions from the Rhone River under certain wind conditions [30], intrusions from the northerly current from the Ligurian Sea [31], and local anthropogenic influences due to the proximity of the Marseille agglomeration (e.g.sewage plants) [32].
In this paper, we aimed to test the hypothesis of a bottom-up relationship between environmental conditions, phytoplankton, and zooplankton in the BoM.To achieve this goal we will first describe a number of zooplankton indices and examine their interannual and seasonal variations in order to identify major changes over the last two decades.We will then relate zooplankton changes to variations in other monitored abiotic and biotic parameters, and discuss how the observed changes might be the result of bottom-up forcing.Finally, we will discuss how this could affect planktivorous fish.

2.1-Environmental data
Sampling was carried out in the Bay of Marseille (N-W Mediterranean Sea) at the Frioul monitoring station (60m depth, 43.24˚N; 5.29˚E) of the long-term national littoral observation program SOMLIT (www.somlit.fr), in the eastern part of the Gulf of Lion (Fig 1).This monitoring survey has generated an extensive collection of physical, chemical, and biological data (see summary of the variables in Table 1).Twice a month since 1994, a vertical CTD-oxygen-fluorometer cast of the whole water column (0-55m) and Niskin bottle sampling were performed at three depths (subsurface -1 m-, bottom -55 m-and fluorescence maximum -variable).The values of each measurement obtained at the three depths were averaged.Climatic and meteorological indices completed the environmental database (Table 1).

2.2-Zooplankton data
In parallel with the SOMLIT monitoring survey, mesozooplankton samples have been collected since 2005 through vertical hauls using a WP2 200 μm mesh size plankton net (water column sampled: 0-55 m).The samples were fractioned in two.The first half of the cod-end content was preserved in a 4% buffered formaldehyde solution for digitalizing; the other half was maintained in cold condition for zooplankton size-fractioned dry weight (biomass) measurements.Complementary information concerning the environmental and zooplankton monitoring is given in the S1 File.
2.2.1-Zooplankton size-fractioned dry weights.The aliquot of cod-end content maintained in marine water was fractioned on a sieve column in 6 size fractions: 80-200 μm, 200-300 μm, 300-500 μm, 500-1000 μm, 1000-2000 μm and >2000 μm.The material from each sieve was recovered on pre-combusted and pre-weighed Whatman GF/F filters and dried for 48h at 60˚C for biomass measurement.Because the WP2 200 μm mesh size plankton net does not efficiently sample the smallest and largest size fractions, the 80-200 μm and >2000 μm fractions were not considered in our study.

2.3-Data analysis
Data analyses were conducted on R software version 4.1.2[40].

2.3.1-Description of zooplankton community structure indices.
Size spectra were constructed using images with ESDs between 200 μm (limit of particle detection by Zooscan) and 1660 μm (i.e.outliers were excluded according to the Tukey method [41]).Following the method described by Nerini and Ghattas [42], the shape of the zooplankton size spectra was studied by means of Functional Principal Component Analysis (FPCA).The FPCA is a statistical tool that ordinates sample units according to the characteristics (i.e. the shape) of their functional entity (i.e.here a curve/size spectrum).This methodology was applied as follows: (i) zooplankton size spectrum densities for all samples were estimated using Kernel density estimator; (ii) the density functions were converted into a continuous functional object, a curve, using 250 B-splines functions; (iii) the curves obtained from all samples were ordinated by means of a FPCA [43].We applied the FPCA on the coefficients of the B-spline expansion model as explained by Pauthenet et al. in [44].These new variables which concentrate the variance of the system are the principal components (PCs) and represent the most significant modes of the data variation.
The sample scores on the first two PCs were used as time series to describe the variations of the size spectrum shapes and by extension the structure of the zooplankton community size.
In addition, multidimensional scaling (MDS) was applied on the species abundance table to disentangle the diversity structure with Manhattan distance.The sample scores of the three MDS axes were used as time series to describe the variations of the diversity structure.
The R packages 'fda' [43] and 'vegan' [45] were respectively used to perform the FPCA and MDS analysis.
2.3.2-Exploringinterannual changes in zooplankton data.Time series data were averaged on a monthly basis and missing values were imputed using the Multiple Imputation with Principal Component Analysis (MIPCA) procedure using the R package 'MissMDA' [46].Trends in monthly time series were highlighted using Local Polynomial Regressions (LOESS, using 75% of the neighborhood data).Breakpoint detection procedure ('strucchange' package in R [47]) was applied to highlight the structural change on every single time-series.Until three breakpoints were investigated, the optimum number of breakpoints was defined according to the Bayesian Information Criterion (BIC).

2.3.3-Zooplankton seasonality.
Level plots were used to describe seasonal patterns in the various zooplankton series: biomass, abundance (see 2.2-Zooplankton data) and structure descriptors (see 3.1-Community structure indices).In order to obtain a smoothed picture of the level plot, LOESS functions were applied over 30 equally spaced points on the year/month  SPECIMED [35] and PHYTOBS [36] Ratio diatoms counts: dinoflagellates counts (subsurface) Ratio Diat:Dino -Environmental 2010-2020 SPECIMED [35] and PHYTOBS [36] High Nucleic Acid Bacteria counts (subsurface) HNA bacteria cells.mL- grids.Seasonal patterns were assessed by fitting sinusoidal models.The seasonal pattern was considered significant when the sinusoidal effect was significant (p-value<0.05).The seasonal onset timing of each zooplankton group was defined as the date when 20% of the yearly cumulative value was reached [8,48].

2.3.4-Investigating relationships between environmental and zooplankton time series.
Dynamic Factor Analysis (DFA, [49]) is a statistical method that aims to find common trends among observational data.In our case DFA was applied in order to find common trends between environmental and zooplankton data.DFA belongs to the family of the statespace models and allows determination of M hidden common trend(s) among N standardized (to reduce the impact of the variable scale on the analysis) data time series, when M<<N.DFA models relate observational time series (y) to the sum of the M hidden trends linear  combination and noise.This can be represented by the following Eq (1): Γ is a matrix of dimensions N x M and contains the factor loadings of the N factors for M detected hidden trends, α t is a vector of length M containing the M common trends at time t.ε t is the noise component which follows a normal distribution of mean 0 and variance-covariance matrix R.
DFA models require user specification for certain parameters: • Firstly, different assumptions on the noise, R, variance-covariance matrix among factors can be made: • same variance and no covariance ('diagonal and equal') • different variances and no covariances ('diagonal and unequal') • same variances and same covariances ('equalvarcov') • different variances and covariances ('unconstrained') • Secondly, the number of M hidden trends must be specified as low as possible to facilitate interpretation (M<<N).
Correction of Akaike Information Criterion (AICc) was applied to select the most likely data model.
Multiple DFA models were computed on: 1. -all environmental data for the 2005-2020 period.From one to three common trends were investigated among the 34 environmental times series with the four different assumptions on the R matrix.A total of 12 models were computed on environmental time series.
2. -zooplankton data for the 2005-2020 period.From one to three common trends were investigated among 11 zooplankton times series with the four different assumptions on the R matrix.A total of 12 models were computed on biological data.
3. -environmental winter conditions (average values between January and March) and dates of zooplankton seasonal onset.From one to three common trends were investigated among the 34 environmental and 19 zooplankton time series with the four different assumptions on the R matrix.A total of 12 models were computed.
As DFA aimed at finding common trends among short time series, in (1) and ( 2), time series were averaged at a biannual time step.This allowed us to reduce the computing time for the models.The series were de-seasonalized (i.e. when seasonal patterns were detected, by means of partial autocorrelation functions) to focus on interannual variations.DFA can be performed with missing data.For models 1 and 3 some time series were shorter (i.e.micro-, pico-, nano-plankton data started in 2009); the results of the models computed for those variables were therefore not discussed before 2009.
Dynamic factor analysis models were computed with the MARSS package [50].

3.1-Community structure indices
The two first PCs of the FPCA performed on the size spectra curves explained 66.7% of the data variability (

3.2-Zooplankton interannual variations and relation with the environmental trend
The zooplankton community in the Bay of Marseille is dominated by calanoids (Fig 3A) which represent more than 50% of the individuals sampled.The second most represented taxa were appendicularians, oithonoids and crustaceans which represent with calanoids, on average, more than 80% of the zooplankton community identified in this study.Obvious changes in relative abundance of some species occurred over the observation period.The dominance of calanoids seemed to have declined from more than 60% before 2012 to 50% after 2017, in contrast to the relative abundance of crustaceans which almost doubled from 5 to 10% during the same period.The changes in the relative abundances were mainly attributed to the diminution of copepod calanoid and oithonoid absolute abundances (Fig 3B).More details concerning the patterns of variation in the zooplankton community structure in species were observable on the interannual series depicted at monthly scale (Fig 4).The abundances of bivalves, chaetognaths, harpacticoids, crustaceans, pteropods, salps, appendicularians, fish eggs, and cnidarians seemed to have increased during the study and breakpoints were detected for all these taxa (except for cnidarians) in the last five years of the time series.The series of the two main copepod taxa (calanoids and oithonoids), copepod nauplii and total zooplankton abundance showed a diminution, and a breakpoint was detected for calanoids and oithonoids series at the end of 2012 (in November and September).

3.3-Zooplankton seasonality and relationship between spring-summer production and winter environment
Most of the zooplankton series presented a clear seasonal pattern (Fig 7) but with different shapes according to groups and size fractions.Seasonality assessments were significant (p-values<0.05)for every variable except for the second axis of FPCA (see S4 File).The total zooplankton abundance and biomass, the four-biomass size-fractions between 200 and 2000 μm, and the sample scores on the first axis of the MDS followed a similar seasonal pattern to calanoids, oithonoids, ergasilida, bivalves, copepod nauplii and harpacticoids with their highest abundance values in spring (April-May).The remaining zooplankton groups (appendicularians, chaetognaths, cnidarians, crustaceans, pteropods, salps, fish eggs) and the two remaining MDS axes seemed to reach their highest values during summer.

4.1-Seasonal and interannual variations of multiple zooplankton indicators
The seasonal patterns of the main zooplankton groups in the Bay of Marseille during the time series are consistent with zooplankton successions observed in Mediterranean coastal areas [5,27,[51][52][53][54], with a zooplankton community dominated by copepods in response to the spring phytoplankton blooms and a second late summer -fall peak of abundance.In the Gulf of Lion, phytoplankton-zooplankton peaks are associated for both periods with strong grazing activity [14].The second peak is associated with higher diversity indices (as shown in Figs 7 and 8) linked to a relative diminution of the copepod contribution (mainly calanoids and oithonoids) and the appearance of other taxonomic groups (either herbivorous taxa feeding on small size phytoplankton, such as salps and pteropods, or carnivorous taxa, such as chaetognaths, crustaceans and cnidarians), a general pattern already described in the NW Mediterranean Sea [11].In addition, this result showed a reciprocal relationship between total abundance/biomass (mostly driven by copepods) and diversity at the seasonal and interannual scales.This relationship was already observed in the Mediterranean Sea and in the global ocean [55,56].
Our results highlight a clear shift in the zooplankton community structure (size and taxonomic assemblage) in 2013.The community structure index, MDS 1 axis, clearly represented the variations in dominance of copepods calanoids and oithonoids compared to other taxa.Since the MDS 1 index aims at summarizing dynamics of multiple taxa, it is obvious that some breakpoints are observed later or earlier when investigating interannual variations of  zooplankton community dynamic at the seasonal and interannual scale (Figs 3 and 4 and 7 and 8).Our results are consistent with the zooplankton dynamics of the Bay of Calvi [27] and the Bay of Villefranche-sur-Mer [54] during the same decade.The main changes in the zooplankton community in the Bay of Marseille observed in 2013 (see Figs 4 and 6) appear to have occurred earlier than similar variations in the whole mesozooplankton community observed after 2015 at Villefranche-sur-Mer [54], where crustaceans (including copepods in their study) showed a diminution of abundance and a lower size spectrum slope (mainly due to the loss of small copepods).In their analysis of multiple time series in the Mediterranean Sea, Berline et al. [5] did not find any common zooplankton interannual patterns at basin scale.Nevertheless, our results and those at Villefranche-sur-mer and Calvi suggest that the coastal regions of the Ligurian Sea might be similarly impacted by the environmental changes and by their connections through the Ligurian current [57].As the re-analysis of time series with recent data proved that some environmental-biological processes relationships may break down [58], reanalyses of zooplankton time series in the Mediterranean Sea with longer time series are required to assess some large-scale processes.
Our results showed different interannual patterns among the size-fraction biomass dynamics (see Fig 4).The sudden decline in 2008 for the largest size-fraction (1000-2000 μm) was not simultaneously observed for the three other fractions (between 200 and 1000 μm) which decreased between 2 and 4 years later.No major change in the community diversity was observed in 2008, while the diminution of the larger size fraction biomass was accompanied by a shift in the size structure of the zooplankton community.Previous studies in the Bay of Marseille have shown that the taxonomic assemblages are not the same in the four size fractions [59,60], impacting differently the biomass dynamics of each size fraction.The size-fractioned approach enabled us to highlight these differences between the size fraction biomass dynamics that would have gone unnoticed if considering the whole biomass.This revealed the strong link between the biomass size-fraction dynamic and community size structure.
The different community level indicators (i.e.size-fraction biomasses, total zooplankton abundance, size-and diversity-structure) used in our study highlighted shifts occurring in different years (Figs 4 and 6) denoting structural changes that would not have been perceived by more aggregated variables (e.g. total abundance or total biomass).In 2013 the assemblage composition from image analysis was sufficient to explain the shift.However, a finer taxonomic reanalysis would certainly be necessary to understand the 2008 shift potentially due to changes within the calanoid group, a group with high diversity in species composition, size distribution and trophic behavior [61,62].In addition, such a more detailed taxonomic identification, up to separation of adult and copepodite stages, could help to better interpret the fine changes in size spectrum, highly sensitive to the underlying assemblage compositions and their associated size distribution.
Overall the multiple indicator approach used in our study helped us to better characterize and understand the zooplankton community dynamic.Pitois et al. [2] also advocated a multiindicator approach as necessary for describing zooplankton structural changes.

4.2-Variations of the environmental conditions and link with zooplankton
In the present study, we aimed to understand the observed zooplankton variations within an ecosystemic framework at the study site (as shown in Fig 5) including changes in the abiotic conditions and in the lower trophic levels (micro-nano-pico-plankton).We described multiple environmental trends: • Higher values of salinity and oxygen between 2005 and 2007 (Fig 5 , trend 1).This anomaly can be interpreted as the signature of stronger offshore deep-water inflows during this period.By analyzing satellite images in those years, Mayot et al. [63] highlighted particularly high phytoplankton productivity induced by strong deep water convection events occurring over a broad area in the Provencal basin with effects extending to surrounding coastal areas, including the Bay of Marseille.
• The reduction in nutrients and particulate matter inputs observed in our time series (Fig 5 ,  trend 1) is consistent with observations already made in the Bay of Marseille, and appears to be a general trend in French Mediterranean coastal waters [64,65].The physical and chemical characteristics of the Bay of Marseille might have been impacted by the diminution of the Rho ˆne nutrient inputs in the Gulf of Lion [66] and/or by the proximity to the Marseille sewage treatment plant where the treatment efficiency has improved since 2008 [32,67].While ammonium, nitrite and phosphate decreased, nitrate concentration was stabilized (and even presented a slight increase in the second part of the study, see S3 File).As autotrophic planktonic groups show preferences in nutrient uptake [68], the variations in the stoichiometry of nutrient sources and their concentrations may have affected their distribution.
Our observations evidenced clear changes in size and abundance of the micro-, nano-, pico-plankton in 2013 which were concomitant with changes in the mesozooplankton community structure (Figs 5B and 6B and 6C).This supports the hypothesis of a regime shift in the whole pelagic food web through bottom-up cascading effects [69,70], similarly to those described in Rietkerk et al. [71].The concomitant decline of calanoid copepods and increase in appendicularians and salps represented an indicator of such change in the zooplankton community.Because appendicularians and salps can feed on smaller particles, the trophic competition with calanoid copepods (feeding on larger particles) seemed to be reduced [72,73].When copepods are unable to feed on large microphytoplankton, they often switch to microzooplankton [74,75] resulting in a higher trophic level of copepods [69].The hypothesis of such a relationship between the phytoplankton size and zooplankton assemblages has been suggested in the NW Mediterranean Sea [54,76].These changes in the Bay of Marseille were accompanied with the appearance of predators such as chaetognaths or cnidarians.In the Catalan Sea, the predation impact of chaetognaths on the standing stock of copepod was low, but chaetognaths can exert a high level of predation on specific copepod species and stages under food-limited conditions [77].A species level taxonomy would help to quantify the top-down pressure of predators on copepods species dynamics.
In parallel, by analyzing the winter environmental conditions, we observed a relationship between NAO, precipitation, temperature, and zooplankton phenology (Fig 9).Although at Mediterranean scale no relationship was found between NAO and zooplankton by Berline et al. [5], several studies suggested that winter conditions at NW Mediterranean Sea regional scale (inducing strong variation in the NAO index) impacted the following zooplankton spring peak [4,27,28], as shown by our results (and discussed in the previous section).Since 2014, winter deep convection (associated with climatic forcing NAO) declined in the NW Mediterranean Sea causing a stratification and water warming [78].This was concomitant with the increase of winter temperature, NAO and precipitation in the Bay of Marseille, and (i) a later seasonal onset of the 200-1000 μm zooplankton biomass, copepods (calanoids and ergasilida), appendicularians, total zooplankton (ii) an earlier seasonal onset of salps.Therefore, we cannot exclude the hypothesis that winter large-scale forcing impacted the coastal zooplankton dynamics in the Bay of Marseille in the 2010s.

4.3-Implications of the variations in the pelagic ecosystem
The conceptual schema, Fig 10, summarizes the interannual changes analyzed in the present paper concerning the environment of the Bay of Marseille and its mesozooplankton community.The results of the Bay of Marseille time series contribute to the documentation of the alterations in the Mediterranean pelagic ecosystem [79], particularly changes occurring in the Gulf of Lion.Espinasse et al. [17], on the basis of a large-scale survey in winter and spring, defined three zooplankton habitats in the Gulf of Lion: the continental shelf, the Rhone influence zone and the Occitan littoral zone.They found that the zooplankton habitat in the BoM presented the same environmental and biological characteristics as most of the Gulf of Lion continental shelf.Therefore, the Bay of Marseille time series can be considered as a sentinel zone for monitoring environmental effects from the open sea as well as the coastal area.
The observed sudden diminution of the larger mesozooplankton size fraction biomass (shown in Fig 4) coincided with the shift in fish body condition at Gulf of Lion scale evidenced in 2008 [20] and with observed changes in the diets of small pelagic fishes [19].This supports the hypothesis of a bottom-up control of the pelagic food chain up to the planktivorous fishes in the Gulf of Lion [24,66], as shown in other areas in the Mediterranean Sea [81,82].
By taking into account indicators related to the zooplankton populations seasonal onset, we show that the process of match/mismatch with small pelagic fish [83] may be a process at work in the Gulf of Lion.An increased mismatch between the spring peak of zooplankton biomass (mainly calanoids) and the growth phase of small pelagic fish [84][85][86] could certainly explain their body condition.In their experimental study, Queiros et al. [87], showed that sardine can display adaptative phenotypic plasticity to food condition changes.Under lower food quantity and quality conditions, the smallest phenotypes experienced lower mortality by starvation than the larger ones during the critical post-reproductive period (i.e. at the end of winter).In this context, a delay in the availability of these preferred fish prey might favor smaller fish phenotypes.In addition, the zooplankton biochemical composition and energy content may vary seasonally [60,88] due to the reorganization of different zooplankton assemblages.Further investigations of the interannual variations of the mesozooplankton quality, both in species composition and biochemical content, are needed to improve our understanding of changes in the trophic environment of the small pelagic fishes.

Fig 1 .
Fig 1. Map of the study site.In the left panel, localization of the Bay of Marseille (red) at basin (NW Mediterranean) and Gulf of Lion scale.In the right panel, localization of the Frioul sampling station in the Bay of Marseille.Color gradient refers to the bathymetry on the left panel.Maps were produced using Natural Earth and France-GeoJSON [33] open access data.https://doi.org/10.1371/journal.pone.0292536.g001 /doi.org/10.1371/journal.pone.0292536.t001 Fig 2A and 2B).In each figure, positive and negative influences of the factorial axes on the size spectra shape are displayed.The variation modes of the size spectra shape on the first axis explained 45.9% of the data variability.Positive values on the first factorial axis (Fig 2A, red curve) were associated with a displacement of size spectra mode towards higher ESD.Negative values on first factorial axis (Fig 2A, blue curve) were associated with a displacement of size spectra mode toward lower values.The second axis of the FPCA explained 20.8% of the data variability.Positive values on the second PC (Fig 2B, red curve) were associated with the displacement of the size spectra peak toward lower ESD and the increase of the proportion of larger zooplankton organisms (>500 μm ESD).Negative values on the second PC (Fig 2B, blue curve) were associated with the displacement of the size spectra peak toward higher ESD and the decrease of the proportion of larger zooplankton organisms (>500 μm

Fig 2 .
Fig 2. Results of the FPCA (A & B) and MDS (C & D) analysis to derive community structure indices.Influence of factorial axis (A for first PC and B for second PC) on the shape of the functional object (size spectra curve).Black line is the mean functional curve of the FPCA.Red curve (+ symbol) and blue curve (-symbol) represent respectively the positive and negative influence of the factorial axis on the size spectra shape.Limits of the x-axis of A & B graphs are set to 200 to 1000 μm, as the main variations of the influence of the factorial axis on the functional object are observable before 1000 μm ESD.Representation of the ordination of the MDS (C: first vs second axis, D: first vs third axis).The black dots correspond to the projections on the 2D space of the samples.The red labels correspond to the centroids position of the taxa.Taxa abbreviations correspond to the following groups: app: appendicularians; biva: bivalves; cala: calanoids copepods; chae: chaetognaths; cnid: cnidarians; crust: crustaceans (without copepods); eggs: fish eggs; erga: ergasilida copepods; harp: harpacticoids copepods; naup: copepods nauplii; oith: oithonids copepods; pter: pteropods; salp: salp.https://doi.org/10.1371/journal.pone.0292536.g002 Abundance estimations of ergasilida presented a breakpoint in April 2007, resulting in an increase in abundance.Interannual variations of the sample scores on the first MDS axis diminished (with a breakpoint detected in April 2014); this observation was based mainly on the diminution of the copepod dominance while abundances of other taxonomic groups such as chaetognaths and salps increased.The second axis highlighted the increase in abundances of bivalves, harpacticoids and pteropods (breakpoint in October 2012).The third axis represented the diminution of nauplii and increase of crustaceans abundance with the breakpoint in April 2017.The smoothing of the first axis of the FPCA revealed that after a displacement toward lower ESD values (between October 2009 -March 2012), the size spectra shifted toward higher ESD after March 2012.No outstanding interannual variations of the second PC were observable and no interannual structural change of the series was detected.The biomass of the 1000-2000 μm size-fraction presented a breakpoint associated with a decrease in September 2008.The three other zooplankton biomass size fractions decreased later, but no breakpoints were detected.On the 12 models performed on the environmental data for the 2005-2020 period that have converged, the DFA model with a diagonal and equal R matrix and displaying two trends was considered as the best data model (see S3 File).The first trend (Fig 5A) was (a) positively characterized by variations in S, concentrations of O 2 , NH 4 , NO 2 , PO 4 , POC, PON; and abundances of bacteria (tot, HNA and LNA), Crypthophycea, Synechoccocus, Prochlorococcus, pico-, nanoeukaryotes and microphytoplankton, and (b) negatively characterized by variations in the size index of Crypthophycea and nanoeukaryotes.This trend decreased from 2005 to 2009, increased suddenly between 2011-2014 and decreased again slowly until 2021.The second trend (Fig 5B) of this model was characterized by a sudden increase between 2013 and 2014.Diatoms: dinoflagellates ratio was positively weighted while size indices of bacteria, Synechococcus, Prochlorococcus and nanoeukaryotes were negatively weighted.The fit of the data model to the environmental time series is represented in S3 File.

Fig 3 .
Fig 3. Year-to-year variations in the A) relative and B) absolute abundance of the 13 best recognized zooplankton taxa.Note that the graph is zoomed on the y values below 0.6 to better see the species with low relative abundance.Abundances in the bottom panel are expressed as number of individuals.m - .The taxa are sorted by importance over the whole zooplankton time series.https://doi.org/10.1371/journal.pone.0292536.g003

Fig 4 .
Fig 4. Representations of the monthly interannual times series of the zooplankton variables (black dots).The red continuous (and dashed) lines correspond to the smoothing (and its 95% confidence interval) by means of local regressions to describe the trends.The blue lines correspond to the mean values of the series before and after a breakpoint date.For series without blue lines, no breakpoints were detected.Series of size-fractions biomasses, in mg.m -3 , (Biom 1000-2000, Biom 500-1000, Biom 300-500, Biom 200-300 μm), taxa abundances and total zooplankton abundance, in individual.m - , are log transformed (log10(x+1)).https://doi.org/10.1371/journal.pone.0292536.g004

Fig 5 .
Fig 5.The two trends depicted by the best DFA model (diagonal and equal R matrix and two trends) on environmental data (2005-2020) and factor weights associated with each trend.Variables with positive weights are associated positively and depicted similar trends to those shown in the left panels.Inversely, variables with negative weights displayed inverse trends.Note that only the most significant factors (i.e. with weight higher than 0.2 in absolute values) are represented.See Table1for abbreviations of factors.https://doi.org/10.1371/journal.pone.0292536.g005

Fig 6 .
Fig 6.The two trends depicted by the best DFA model (diagonal and unequal R matrix) on zooplankton data (biomass, size and diversity structure indices and total zooplankton abundance series) between 2005-2020 and factor weights associated with each trend.Variables with positive weights are associated positively and depicted similar trends to those shown in the left panels.Inversely, variables with negative weights displayed inverse trends.A) and B) represent respectively the first and second zooplankton trends with their factor weights.C) trend 2 of zooplankton (in red) and trend 2 of environment (in black, see Fig 5B) are represented together.Vertical dotted blue line corresponds to the common breakpoint detected between both trends.https://doi.org/10.1371/journal.pone.0292536.g006

Fig 8 .Fig 9 .
Fig 8. Level plots of the monthly interannual times series of the zooplankton variables.Circle colors correspond to the month value and the background color corresponds to smoothed values with LOESS function over 30 equi-spaced points.Series of biomass size-fractions, in mg of dry weight.m - , (Biom 1000-2000, Biom 500-1000, Biom 300-500, Biom 200-300 μm, Total Biomass), taxa abundances and total zooplankton abundance, in individuals.m - , are log transformed (log10(x+1)).Yearto-year variations of the date of seasonal onset, considered as 20% of annually cumulative value for biomass and abundance time series, are represented by the blue lines.https://doi.org/10.1371/journal.pone.0292536.g008

Fig 10 .
Fig 10.Summary of our main results by means of the revisited Ramon Margalef mandala [80].In the schema, the significant interannual changes of the zooplankton community and environmental parameters are represented along the timeline materialized by the black arrow.The variables around the green square exerted a direct or indirect bottom-up control on the interannual zooplankton community dynamic: (1) Nutrients stoichiometry, (2) POM and nutrients quantity, (3) Micro-, Nano-, Pico-plankton size, (4) Winter conditions.The taxonomic groups within in the central square highlight the main variations of the zooplankton community composition.The variations in the zooplankton community traits are represented around the red square: (5) Biomass and abundance, (6) Size structure and diversity, (7) Seasonal onset.Finally, the implications for planktivorous fish are represented within the blue squares: (8) Quantity of food variation, (9) Preferred prey variation, (10) Match-Mismatch.https://doi.org/10.1371/journal.pone.0292536.g010